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Abstract 

In this paper is presented the new modelling framework combining the repli- 
cator dynamics (which is the standard model of the frequency dependent selec- 
tion) with the Leslie Matrix model of the age-structured population. Firstly the 
continuous version of the discrete Leslie Matrix model is derived. It is shown 
that Euler-Lotka equation is satisfied when new model reaches the steady state 
(ie. stable related frequencies between the age classes). Due to the long expected 
lifespan of an individual in comparison with the ecological timescale, the real 
life model should contain a large number of equations. This problem is solved 
by the introduction of the large age classes concept. The underlying assumption 
is that within a single large age class the individuals do not differ in the de- 
mographic parameters (fertility and mortality). Then according to this result, 
a more complex model containing different individual strategies is presented. 
The methodology of the multipopulation games is used for derivation of two, 
mutually equivalent systems of equations. First contains equations describing 
the evolution of the strategy frequencies in the whole population completed by 
subsystems of equations describing the evolution of the age structure for each 
strategy. Second system contains equations describing the changes of general 
populations age structure, completed with subsystems of equations describing 
the selection of the strategies within each age class. 

Introduction 

Among the most important approaches to modeling of evolutionary pro- 
cesses are life history optimization and evolutionary game theory. The classical 
approach to life history theory (Stearns 1992, Roff 1992) relies on optimiza- 
tion models. In effect, there are no interactions among individuals or density 
dependence: 

"Life history evolution usually ignores density and frequency dependence. 
The justification is convenience, not logic, or realism" (Stearns 1992)" 

On the other hand, in classical game theoretic models there is no age or 
stage structure. Payoffs describes the averaged lifetime activity of an individual, 
which can be found for example in Cressman 1992: 

"...an individual's strategy is fixed over lifetime or, alternatively, the life 
history of an individual is its strategy." . 

The classical approach to evolutionary game theory that can be found in the 
textbooks (Maynard Smith 1982, Cressman 1992, Hofbaucr and Sigmund 1988 
and 1998, Weibull 1995, Nowak 2006) mostly relies on the mathematical struc- 
ture of a two-person matrix game. However, during recent years there has been 
a huge development in a field, which extends the matrix game paradigm de- 
scribing population outcomes of local two-player interactions. There are papers 
related to modeling of hierarchy formation (Broom (2002), Broom & Cannings 
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(2002)) where more than two individuals are involved in a single interaction 
event and density dependent selection (Krivan et al. 2008, Argasinski and 
Kozlowski 2008). In addition, there are generalizations of basic evolutionary 
game approaches to cases where individuals are characterized by a continuous 
trait or a set of continuous traits. Therefore, there is progress in game theoretic 
methods to using similar methodology to that of life history theory. The clos- 
est approaches are state based games (Houston and McNamara , 1988, 1991, 
1992, 1999, Webb et al, 1999, McNamara et al., 1994, McNamara et al, 2004, 
McNamara et al., 2009). Another important approach in theoretical evolution- 
ary ecology is individual based modeling (multiagent computer simulation tech- 
niques). This is a methodology based on the assumption that the individual is 
the basic subject of selection and that a model should explicitly consider unique 
traits of individuals. The creators of this approach (Lomnicki 1988, Grimm, 
Wyszomirski, Aikman, Uchmanski 2003, Grimm and Railsback, 2005) criticize 
classical mathematical models as oversimplified. 

Methods used in life history optimization are closely related to classical de- 
mographic methods and are expressed in demographic terms. However, how to 
construct a general description of the relationships between demographic struc- 
ture and population dynamics is still unsolved problem (Caswell 2011). Game 
theoretic models operate in abstract terms of costs and benefits measured in 
units of fitness mostly without deeper insight in their meaning or interpretation. 
This problem was analyzed in Argasinski and Broom (2012) where relationships 
between classical demography and evolutionary games are described in detail. 
The main conclusion there is that instead of the excess from average fitness, 
models should be described explicitly by mortality and fertility, which are ba- 
sic opposite forces shaping population dynamics. Those results are sigificant 
progress in ecological realism; they emphasize the role of background mortality 
and fertility or turnover of individuals (this was originated by Argasinski and 
Kozlowski, 2008). However, that approach is still very primitive. Background 
mortality is described by a single constant, which describes an exponential de- 
cay of the population, which implies that the length of an individuals lifetime 
is potentially unbounded, and there is no aging and no age specific payoffs. 
The goal of this paper is to fill this gap and develop a mathematical structure 
combining selection of individual strategies with an age structured population. 

0.1 Classical theory 

The classical approach to the modelling of age structured populations is related 
to Leslie matrices, followin g the following matrix equation: 
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This age-structured growth model suggests a steady-state, or stable, age- 
structure and growth rate. The growth rate can be calculated from characteristic 
polynomial of the Leslie Matrix called the Euler-Lotka equation: 

m i — 1 

fo + £ e~"7, II Sj = 1. 

where r is the intrinsic growth rate of the population. Now let us focus on 
populations diversified by individual strategies. 

When a population with a finite number of individual strategies is consid- 
ered, then a system of differential equations called the replicator dynamics can 
be defined. It describes changes of the population state in time and can be 
derived in the following way. Assume that we have a finite number / of arbi- 
trarily chosen strategies. For each strategy some payoff function r, is assigned 
(for example in matrix form, as in the classical Hawk-Dove Game, however the 
form of payoff function depends on the modelled problem and can be more com- 
plicated). Then the growth of the population of i-strategists can be described 
by the Malthusian equation 

hi = run. (1) 
Then by following a change of coordinates, 

ft = - (2) 
n 

where n = Ei=i n i 1S the population size, we can derive a system of ordinary 
differential equations (see e.g. Cressman (1992)), 

qi = qi(n-f) (3) 

where f = Ei=i ft r « * s the average payoff in the population. However, instead of 
the Malthusian parameter describing payoff we can explicitly consider individual 
fertility fa and individual mortality di of an i-strategist. In effect we obtain the 
following variant of the replicator equations that can be called the sex and 
violence equations, 

4< = Qiifi -f-ck + d) = qi ((fi - f) - (di - d)), (4) 

for details see Argasinski and Broom (2012). It was shown (Argasinski, 
2006) that every single population system described by the replicator equations 
(|3"T) can be divided into the product of subsystems describing the dynamics in 
arbitrary chosen subpopulations (described by a vector of frequencies q l for the i- 
th subpopulation) and an additional system describing dynamics of proportions 
between those subpopulations p (see appendix A for details) . The dynamics in 
each subpopulation will have the form ((3]) and will depend on the excess of the 
strategy payoff from the average payoff in this subpopulation. Therefore, the 
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same operation can be carried out for equations (U]). In effect we obtain the 
following system: 



'/;((/: n k <?))> 



(5) 




in the i- th subpopulation and f l and d % are the mean fertility and mortality, 
respectively, in the i- th subpopulation and / and d are the respective values 
in the global population. When fertility and mortality are density dependent, 
the system should be completed by an equation on the population size (scaling 
parameter) : 



0.2 A continuous version of the Leslie Matrix 

A major technical difference between Leslie matrices and replicator dynamics 
is the fact that the first describes discrete system and the second is a system 
of continuous differential equations. The Leslie matrix can be presented as a 
system of algebraic equations, describing changes of age structure with m age 
classess: 

n (t + 1) = f n (t) + /irc-i (*)••• + f m n m {t), 
rii(t + 1) = Sj_ini_i(t), 
n m (t + I) = s m -xn m -x(t). 

We can imagine a continuous version of this system: 



To solve this system we can apply the theory of delayed differential equations 
because equation (J5|) leads to rii(t) — no(t — iT)sQS\S2...Si-i where r is the 
respective time interval. This will lead to Euler-Lotka equation for continuous 
system. However we will use alternative way to do it by application of the 
replicator dynamics. What happens when we describe the system in related 



h = n(f - d). 



(7) 



hi = s n - m, 



(8) 
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frequencies Oj = n i/^2j n j an d a scaling parameter (see Appendix b)? We 
obtain the following system: 



di-iSi^i - ai \ a * (/» + s i)j for i = !) -i m > ( 9 ) 
£oi(/i + »i)-l , (10) 



\i=0 



where ao = 1— and the malthusian parameter is r = X)i=i a « (/•» + s i)~ 

1. For a detailed derivation see Appendix B. For the system ([9]) the following 
steady state can be calculated: 



a s 

a\ = 



a _ CLqSqSi 

(EiLo a i(fi + S i)) 



i-1 



(11) 
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3=0 

ai 



(E™o a *(^ + S 0) 1 



In addition in state [TT] Euler-Lotka eqution is satisfied (for a derivation and 
proof, see Appendix C) The above model is the continuous equivalent of the 
basic approach related to the Leslie matrices used in life history modelling. 
According to this theory natural selection promote strategies which maximize r 
(which is the root of the Euler -Lotka equation). However, in the above model 
there is no explicit density dependence. It is hard to imagine the unlimited 
growth of the population indicated by this model. In addition, the problem 
of a universal fitness measure is more complicated when density dependece is 
explicitly considered. Thus our model should be develop further. 



0.3 The impact of density dependence 

To explicitly model density dependence we can apply a selectively neutral mech- 
anistic approach from evolutionary games (Argasinski and Kozlowski, 2008; Ar- 
gasinski and Broom, 2012) which is related to a modification of the classical 
logistic growth (Kozlowski 1980, Kuno 1991, Ginzburg 1992, Gabriel 2005, Hui 
2006, Argasinski and Kozlowski 2008). This constitutes the eco-evolutionary 
feedback affecting the fitness (Argasinski and Kozlowski, 2008; Zhang and Hui, 
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2011; Argasinski and Broom, 2012; Argasinski and Broom, 2013b). This ap- 
proach is closer to current developments in the evolutionary theory focused on 
the relationship of selection processes with ecological factors (Schoener 2011, 
Morris 2011, Pelletier.et al. 2009). In this approach each fecundity should be 

multiplied by the logistic supression coefficient ^1 — —^J describing the survival 

of newborns introduced to the environment (where the carrying capacity K is 
interpreted as the maximal environmental load, Hui 2006). As was shown in 
Argasinski and Broom (2013), at population size equilibrium density depen- 
dent growth equation can be expressed in terms of the turnover coefficient (i.e. 
the number of newborns per single dead individual during At). If there are 
other strategies in the population, this induces a frequency dependent selection 
mechanism called the nest site lottery in Argasinski and Broom (2013). 

The same reasoning can be applied to the age structured model. The above 
model describes a homogenous population divided only by age classes. Thus 
according to frequency dependent selection, a corresponding model of this type 
should be presented for each strategy present in the population, since frequency 
dependence needs information about the state of the whole population. How- 
ever, in the case of the nest site lottery mechanism, the value of the turnover 
coefficient for the general population is sufficient. Thus our model will be ex- 
tended to the density dependent case by the assumption that there are other 
strategies in the population (but not explicitly considered in the model) and 
their presence affects the dynamics of our model via a phenomenological pa- 
rameter describing average turnover coefficient in the general population. 

Now, let us calculate the turnover coefficient describing the number of new- 
borns introduced into the population per single dead individual removed from 
the population during At. Then F = 2i=o a ifi ^ s the total number of newborns 
with the focal strategy introduced to the population and T = X)"=o a « — Si ) 
is the number of dead individuals of the focal strategy removed from the pop- 
ulation (Argasinski and Broom, 2013). Assume that the subpopulation of the 
carriers of the focal strategy is the part of some general population described 
phenomenologically by total population is of size n pop and the average turnover 
F 

coefficient is L pop = pop (where F pop is the total number of newborns and T pop 

pop 

is the total number of dead individuals in the whole population during At). 
Thus the turnover coefficient for some individual strategy s will be: 

L {f{t), *(*)) = £ = ^JP^f 1 . (12) 

and the equation on the scaling parameter for equation (I10[) describing the 
growth of the population of carriers of this strategy can be presented as 
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/ m rn \ 

" = « (l-^E^-E ^ 1 -*) (13) 

\ i=o i=0 / 

- "S-< 1 ->( -'J (14) 

m 

= „(L(/(t), a (t))(l-to)- l)]>>(l- Si ) (15) 

i=0 

However, here we have an age structured population. When we assume that 
the population is at the stable age structure, then the turnover coefficient has 
the form: 

L (/(*), s(t)) = (16) 

Ei= n ^(l+rr-^i-^) 

3=0 

and the growth equation can be presented as follows (see Appendix D for a 
detailed derivation): 



m i— 1 

n = na (l (/(*),«(«)) (l - ^) - l) ^ J] ^* (1 - „) . (17) 

i=l j=0 U + r ) 

It was shown in (Argasinski and Broom 2013) that when the population 
reaches a stable size then frequency dependent selection emerges on individ- 
ual strategies. This is the nest site lottery mechanism (Argasinski and Broom 
2013) related to the turnover of individuals. Each newborn introduced to the 
environment should find a nesting place. At the population equilibrium size all 
newborns produced by all strategies compete for free places released by dead 
adult individuals. In effect each newborn can be drawn from the pool of candi- 
dates with probability Fj F pop to replace a dead adult and during time At there 
will be T trials. Thus survival probability of a single newborn equals the num- 
ber of all dead individuals during time At (in the number of free places) divided 
by number of newborns introduced (number of candidates) i.e. which becomes 

the reciprocal of the turnover coefficient of the whole population . Thus 

we have: 

n = n °[ — f — m — Ell 77— -Tit 1 -*)- ( 18 ) 

V L P op{t) J i=1 j=Q (1 + r) 

Thus for an L-maximizing strategy it is profitable to maximize T = E»=i Jl (1 — Sj). 

3=0 (1 + rf 
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i-1 

Note that for r = L reduces to the lifetime reproduction Y^T=o II s jfi' 

3=0 

i 

because YlT=o II s j = 1- The structure from (Argasinski and Broom 2013) 

3=0 

works and juvenile mortality will be 1/L (^f(t),s(t)j. However, at (Argasinski 

and Broom 2013) strategies are constant two-dimensional vectors, while here 
we have age depended fecundity and survival. Thus the turnover coefficient 
depends on the age structure of the carrier population of the respective strategy. 
Arguments of L(f(t),s(t)) are sums or integrals £™ ai /« an d a ' s ' 
the nest site lottery works in the same way. Thus at first stage the turnover 
coefficient 

l (/(t), *(*)) -r^ (19) 

EZoUsjil + rr-Hl-Si) 

3=0 

should be maximized. Subsequently, among L-maximizers 

3=0 

should be maximized. Note that younger age classes gave greater weight 
(1 + r) m 1 for L and (1 + r) 1 for T. The average lifetime is 

i-1 

£™i II sj (1 - Si)i, 

3=0 

and so increase in T will lead to decrease in the average lifetime. 

To use the nest site lottery mechanism, we need information about the 
turnover coefficient of the whole population L pop , in this section assumed phc- 
nomenologically. To provide this information we need the model of the age 
structure of each strategy in the population. In effect the model will be ex- 
tremely complicated. However, even the age structured model of the homoge- 
nous population is too complicated to be practical. Every time moment At in 
the lifetime of the individual should be described by a separate equation. A 
realistic model would contain hundreds or thousands of equations. This should 
be simplified before diversity among strategies is introduced into our framework. 
This is the task for the next section. 



1 The behavioural and demographic timescales 

In game theoretic models there is a distinction between behavioural and evo- 
lutionary timescales (Vincent and Brown 2005, Krivan and Cressman 2009). 
There is the possible situation when the behavioural process operates on a differ- 
ent timescale to life history (demographic) dynamics. Generally, we can assume 
that demographic dynamics is shaped by aggregated outcomes of interaction 
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events that occur on the behavioural timescale. For example during a few days 
individual interactions may cause the population to significantly change, while 
differences between different age classes may be measured in years. During the 
time period referred to as an age class, an individual can perform several mat- 
ing events and other activities with a mortality risk. For technical reasons we 
shall define "age class" as a long time period divided into small age sub-classes 
operating in the behavioural time-scale when an interaction may occur. Assume 
that probabilities of death and reproduction are constant during long age class, 
which are sub-divided into small age sub-classes where individuals do not differ 
in payoffs. Under this assumption we can significantly reduce the number of 
equations by considering only the large age classes, with N v = n g + ... + nk, 
where small age class g is the first, and small age class k is the last, of the 
sub-classes of larger age class v. We can assume that all small age classes inside 
long age class are in equilibrium (fTTj) . In effect system © has the form (for a 
detailed derivation see Appendix E): 



(20) 



No = N ^2^fi-(l-s )-D( S0 ^j, 
N v = N v (^j^D (sv-i) - (1 - sv) - D (s v ) ) . (21.) 



where 
D (s v ) 



i=l 

describes the fraction of individuals passing to the next large age class (h v is 
the number of small age classes within the v-th long age class, i.e. its "length"). 



1.1 Frequency dynamics for large age classes 

Now we can transform system (f2Tj) to the related frequencies A$ = . Obtained 

system will be analogous to (|9I10|) . but described in coordinates Ai instead of 
di. All parameters of the system, such growth r should be described in the 
new coordinates. Thus, when we average malthusian coefficients from equations 
(I20l2ip we obtain r = Y^iLi (fi + Si) — D(s m ) — 1, which is equivalent to the 
r from the previous section (D(s m ) describes dying survivors of the last m-th 
age class). In effect we obtain: 
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Ai — Ai 



A,. 



D ( Si _i) - (1 - Si) - D ( Si ) -J^Ai {fi + Si ) + D{s m ) + 1 [2a) 



Ai-iD (si-i) +Ai 



Si - D ( Si ) - ^ {fi + Si ) + D(s m ) 



(23) 



Analogously to equations ([9]) we have A = YliLi Ai and equation on scaling 
parameter: 



N = N I J^Ai {fi + Si ) - D{s m ) - 1 



?=i 



Equilibrium of this system should satisfy the following condition: 



(24) 



Ai 



A i - 1 D{si- 1 ) 



A i ^ 1 D{s i -i) 



D { Si ) + d t +r D { Si ) - Si + YZi A i (fi + s i) ~ D M 
which can be presented in the form: 



(25) 



A 



AoUD {s z ) 

z=0 



II [D (s z ) - s z + i A z (/, + s z ) + D{s m )] 



Since the sum of the Ai terms equals 1, we have: 



(26) 



nD{ Sz 



i + E™ i 



n [d{s z )-s 2 

z=l 



E7=i Az {fz + s z ) + D{s m )] 



(27) 



1.2 Extension to multipopulation replicator dynamics 

Now we can incorporate the above model into a multipopulation evolutionary 
game (see Appendix A for details of the applied methods). Assune that we 
have w strategies and m + 1 age classes indexed from to m. Assume that the 
upper index describes the number of age class and the lower index describes 
the number of the strategy, p describes the strategy (phenotype) fraction and a 
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describes the related size of the age class. Thus, /j and describe, respectively, 
the fertility and survival of the j-strategist in age class i. Two perspectives are 
possible (see fig 1): 

a) First problem is the impact of the age structure in subppopulations 
strategically homogenous on selection of the strategies. It can be described 
by coordinates: 

A®,..., A™ for x = l,...,m the age structure of the ^-strategists 
subpopulation, 

Pi, ....,p w the strategy frequencies in the whole population 

b) The second approach is focused on how strategic selection within each 
age class affects the age structure of the overall population. It can be described 
by coordinates: 

pf , ■■■■,Pw for x = 1, ...,m strategy frequencies in age class x, 
A 1 ,..., A m the age structure of the population 

Now we describe the transition of coordinates between both formulations. 
First we should define the auxiliary canonical metasimplex coordinates (see 
appendix A which is a single simplex without division into subclasses): 

qf = A lP * =PiAf. 

Now according to the Appendix A we can define transitions between the two 
systems: 
a to b: 



Afpi 



A^p-u 



V A x n ' V A x n 
2^2=1 A zPz 2-iz=\ A zP 



A = 



a^pz, a" 



2=1 



z=l 



(28) 
(29) 



and b to a: 



[A x , A" 1 ] 



a°p° x 



A m p™ 



P = 



V™ A z n z ' "' V m A z n z 
l^z=a 71 Px 2^2=0 11 Px 

m m 

J2a z pI,...,^a z p; 



(30) 
(31) 



.2=0 2 = 

Now let us derive systems of equations operating in both coordinates. 
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For system a (see Appendix F for detailed derivation): 



A) = 4-^(4- i ) + 4[4- J D(4)-^(i-^) + « j - J D( S ™ ) )j^) 

Pj = Pi((/ J --/)(l-^)+(5 J --5)-(^7' £, (*r)-^o5a I )). °r(33) 
Pi = Pi((/-/l(l-^)-(^-^-(^( s ™)-^o6a;)), (34) 
and the equation on the scaling parameter 

N = N (7 (l-^j+s- D™ obal lj , (35) 
where ft = Y7=i A) ft, J =Y%=iPifi, *i =YZi*i*)> 

s = Eti Pisu D™ obal = EJ =1 PjAfD(sf). 
For system b (see Appendix G for a detailed derivation): 



\ i i v / ^ ' 



(36) 



^[(-S-s )-^-?)- 50 )] 
4»- 

in 4 - 1 4-4* .^-D'-ls-i-fli- 



pj = ^(p}- 1 ^(4" 1 )-^- 1 )+^[(4-^)- (^( S j)-^)].(37) 



A 4 = A D + A 1 



(38) 



and the equation on the scaling parameter 



where 

^ = E"=iP}^ (sf 1 ) & = E^iPlD (si) D™ obal = A m D r ' 



s * — ~l2v=i Pv s v s — ^2 v =i A v s" 

r = E:=iPitt /=E™i^r- 

Note that D % is the average number of survivors of the i-th age class per 
single i aged individual, while D^ oba( is the average number of individuals that 
survived the last, ra-th age class (and will die afterwards) per single individual in 
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the whole population. Both systems have the same equation on scaling param- 
eter which is equivalent to (|24[) . In effect in both cases equilibrium population 
size is: 



N 



1 + D 



global 



— S 



1 



K 



f 



N = K 1 



1 + D 



global 



f 



— s 



) 



1.3 Background fitness and neutral density dependence 

In this section let us consider the selective neutral (i.e. the same for all strategies 
) factors and their impact on the dynamics. One of such factors is the neutral 
density dependence (Argasinski Kozlowski 2008, Argasinski Broom 2012, Ar- 
gasinski Broom 2013a, Argasinski Broom 2013b). According to this approach 
we can assume that each newborn should find a nest site to live and reproduce. 
The probability of finding a nest site (in effect newborn survival probability) is 
proportional to the fraction of free places and can be described by logistic su- 



pression coefficient 1 — — I , where K is the carrying capacity (number of nest 



places) . Thus each newborn can survive with probability 1 ) . Therefore 



every fecundity should be multiplied by 1 — ) . The case of background fer- 



tility and mortality is more complicated. In (Argasinski Broom 2013b) a more 
realistic approach to the background fitness (fertility and mortality caused by 
other events than the game interactions). It relies on the assumption that some 
number of newborns can be produced and some number of adults can be killed 
during other activities (i.e. background events). We can assume that there is 
a background fertility and mortality acting at the same rate for all individuals, 
in this case those factors will appear only in the equation on scaling parameter. 
Assumption that there are factors acting in the same way on the individual in 
the every age is biologically unrealistic. It is more realistic to assume that out- 
comes of the background events should be the same for all strategies, but they 
should be different for different age classes. Thus we can assume that within 



each age class there is background fertility $ l ( 1 I and background mor- 



tality ^ l . According to the equations (|4l5l6p each additive factor included to 
the growth equation ([1]) , in replicator equations (J3j> will appear as the bracketed 
term describing the excess of the value of this factor from the averaged value 
of that factor in the population. Respectively, after division of the population 
into subpopulations, the intraspecific dynamics will contain the bracketed term 
describing the the excess of the factor value r from the averaged factor valuer 
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in that subpopulation (see [SJ , while dynamics of the subpopulation relative size 
will contain term describing the excess of the avereage factor value in the re- 
spective subpopulation from the average factor value in whole population ([5]). 
Thus according to these relations we can derive the terms describing the impact 
of the factor $' — ty 1 ,i.e. age specific background fertility and mortality on the 
dynamics. 

Let us begin from the system a: 

/ N- 

N = N [J2 ^fi ~ (1 - so) - D (so) 

Nv = N v (j^D («„_!) - (1 - s v ) - D (s v ) 

In subsystem (13"!?|) there are no equation on the age class A® which is influ- 
enced by aggregated fertility Y^iLi thus additional term will be 

m / m m \ 

- A \ ( $i - * l ) = - E - E A )^ 1 ( 40 ) 



In subsystem ([33 



J2 A )&-Y, A W-Y,p>- E4^-E4*M ( 41 ) 



2=1 / 

w m 



:1 j 1=1 I \i=l j=l 1 = 1 



and in equation on scaling parameter (1351) 

/ m tci \ 

Eft E4 $i -E4*M ( 43 ) 



In the case of system b. 

For equations (|36l) individual background fitness can be described as average 

Nj 

per capita number of introduced newborns — ^ Y^—i Aj$ l 

3 

per single newborn present in the age class 0, and specific newborn back- 
TV 1 XT™ A z p z 
ground mortality Vf . Since according to (|30|) — ^ = — ^ = — z ~ Q — - and per 

N j A j A Pj 

capita number of newborns will be: 
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N - 

vffE^*' ( 44 ) 



3 

the bracketed term for equations (|56")l will be 



m A i v i 

§^t** <45) 



771 Ai i w / m Ai i \ 

E>-«°-E«4 E^-*° (46) 

i=l ^ fc=l \i=l ^ fe / 

- E^'-EE^ 1 (47) 

i=l F J fc=l i=l 

Respectively for equations (|37|l there are no direct impact of fertility and we 
have specific age dependent mortality vp 1 the same for all strategies, thus it is 
selectively neutral (the bracketed term is 0). For equations (|38|) we have impact 
of aggregated fertility only in newborns age class and in the next age classes 
acts only specific age dependent mortality "JJ*. Since class is not explicitly 
described in the subsystem (|38|) . the bracketed term will be: 



* l - ^° Erf (E - *°) - E ) (48) 

(w m m \ 

EE^p^-E' 1 '*' (4 9 ) 
fc=l i=l i=0 / 

Cm m \ 

E^^-E^** (50) 
i=l i=0 / 

In equation on scaling parameter ()39|) should be added: 

m m 

E A l <& 1 - ^ A*** (51) 

i=l i=0 

To consider the neutral density dependence fertility factors $ l should be 



multiplied by logistic coefficient ^1 ^ . 

But in the equations we have also operator D. Another decay chanel driven 
by background mortality pressure should be represented in the argument of 
D. Thus we should calculate the per capita fraction of survivors introduced 
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to the next short age class. The problem is that in this case we have two 
types of mortality operating at different rates, thus we cannot reduce this to 
the simple substraction from 1. Let us consider a simplified equation containing 
only mortality components h = n [t± (1 — s) +Ti(l — b))] 

Thus we can use first order Taylor expansion to derive approximation of the 
fraction of survivors passing to the next age class: 

n(t + At) = n(t ) + n(t )At = n(t ) - n(t ) [n (1 - s) +r 1 (l - bj)} At 



n(t ) 1 - 



+ —(1-6) 

Tl 



7"! At 



When TiAt is set to 1 by change of the timescale t = t/ViZ, then we obtain: 



n(t a + At) = n(t ) 1 



1 -*<) + — (I-62) 

Tl 



T 2 



6 2 ) 



And fraction of survivors can be approximated by Sj — \f ,t where VP = — (1 

Tl 

After insertion of the terms P0I BS1 H31 system a will have form: 
system a: 



At = A i ~ 1 D (sT 1 



* l ) + A) [s] - ^ - D (s) - 



(/;-/)+ E4 $i -E^E4* 
\i=i j i=i 

- s) - (A^£>( S ™ - _ D™ lobal ) 



N 
K 



w m 



E^-E^E^ 



i=l 



and equation on scaling parameter 



A\ 



A 
A 



a} - AfDis™ - ^ m ) - 1 



^ / 1 



A 
A 



A 



(52) 



(53) 
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system b: 

For system b (see appendix F for detailed derivation): 



1 

A 5 



m w 



E A % fi - p°i E A * E M )+ E At p) v - p°j E Ai E ri*' 



v.i=l 



i=l fc=l 



iV 



p = 



A 1 



A i ~ 1 D i ~ 1 



~D l - ^s + f(l-^j -A m D m yj - ^E A<$< ( 1_ 



AT 
1? 



bo indeksowane od 1 to nie ma zerowej klasy i fertility per capita 
and equation on scaling parameter 

n = n(j+s- D™ obal i + Y™ i ^ (i - £) - E™ ^) 



1.3.1 Discussion 



In this work was presented a new modelling framework combining evolutionary 
dynamics with demographic structure. This approach can be useful tool in 
the development of the synthesis between evolutionary game theory and life 
history theory. The obtained results clearly show that life cycle perspective 
plays crucial role in the evolutionary processes. In the classical approaches to 
the evolutionary game theory individuals cannot change the properties during 
the lifetime. Thus their life history is the memoryless process, survival of the 
single interaction does not change the state of the individual. This is caused by 
the fact that the classical approaches to the evolutionary games are focused on 
the strategies interpreted as the patterns of behaviour, not on individuals itself. 
The exception from this rule is the state based approach (Houston McNamara 
1999). The new approach presented in this paper is the natural background for 
the dynamic extension of the state based approach. 
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Appendix A Change of coordinates in the space of population states. 



Assume that we want to decompose an entire population into z subgroups. 
Define 

d* = $,...,<] 

as a vector of indexes of strategies exhibited by individuals from i-th sub- 
group (dj S {l,...,u}, Ui number of strategies in the i-th subgroup). For exam- 
ple notation d 2 = [1,3,5] means, that in second subgroup there are individuals 
with strategies 1, 3 and 5. Every strategy should belong to a single unique sub- 
group (and cannot belong to two). Then according to Argasinski (2006), by the 
following change of coordinates: 



(Id' 



for i — 1, z 



we obtain distribution of relative frequencies of strategies in the i-th sub- 
population. Distribution of proportions between subpopulations has the form: 



P = [Pi,--,Pz] = E"=i?4=-)Er=i%? 



where pi is the proportion of the i-th subpopulation. Every decomposition 
into subpopulations can be reduced again to a single population model by the 
opposite change of coordinates q{p,q 1 , ■ ■■ 1 q z ) where: 



Qd* = Pilj 

Note that we can decompose an entire population to z subpopulations. 
When we apply the above transformations to replicator equations, we obtain 
a set of equations that describes dynamics inside subpopulations (intraspecific 
dynamics). Because the argument of a fitness function is a set of relative fre- 
quencies of all individuals q (without division to subpopulations), the opposite 
change of coordinates q(p, q 1 , q z ) should be applied (Argasinski 2006). When 
the set of strategies in each subpopulation is characterized by vector of indexes 
d l , then the system of replicator equations will be: 



q) = q)[W}(q{p,q\...,q z ))-W i {q[p,q\...,q z ))} (57) 

for j = 1, . .., m — 1 and i = 1 , . . . , z 

p s = p s [W s (q{p,q\...,q z ))-W{q{p,q\...,q z ))\ 

for s = l,...,z- 1 (58) 

where W s (q) = Yn=i iV^i '{liP^ 1 > ■■■tQ z )) is mean fitness in the s-th sub- 
population. In simplified notation this system has form (I5I6P In practical appli- 
cations of this method to modeling of biological problems, replicator equations 
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can be defined on decomposed population. This will simplify the formulation of 
the model, because when strategies are initially assigned to subpopulations, then 
there is no need to change their indexes. Choice of subpopulations is arbitrary 
and depends on the biological assumptions underlying the analyzed problem. 
The entire population may be divided into two competing subpopulations of 
hosts and parasites or preys and predators. On the other hand, it may be di- 
vided into two subpopulations of males and females, then interspecific dynamics 
will describe the evolution of secondary sex ratio, and intraspecific dynamics will 
describe changes of frequencies of strategies inside male and female subpopu- 
lation. The entire population can be divided into more subpopulations than 
two. The subpopulations can be divided into subsubpopulations, and the en- 
tire population may be transformed into a complex multilevel cluster structure. 
However, all these structures are equivalent to a single population replicator 
dynamics model. 



Appendix B Derivation of frequency equations 

Initial system [5] can be presented in the form of malthusian equations: 

Ep ( x~\fn n ifi 1 
riiU - n = n }^ i=0 1 

V no 

( n s 

ri\ — sqUq — ni = n\ 1 

V n i 

_ / ^m-l^m-l 



therefore Malthusian growth coefficients Mi — m 1 — 1 — 1 of those equa- 
tions can be presented as frequency dependent functions of age structure, where 
cii = — where (n = are frequencies among age classes. In effect we 

obtain: 

E i= o - — 1 

di-lSi-l 

rti = Hi | 1 

"m-l°m-l 



(I 



n i 



therefore, this system can be presented as system of frequency dependent 
replicator equations hi = ai (Mi — r) and a single equation on scaling parameter 
n = nr. what we need is an average malthusian growth rate: 
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(V^ra a iji -, \ . v^ m I Ci—iSi — 1 \ 

= 127=0 a *fi - a o + YT=i K-iSi-i - «*) 

= J2T=0 a ifi+YT=l Oi-lSi-l-l = E™0 a » (/* + S i) _1 

therefore we can formulate system of frequency dependent replicator equa- 
tions by transforming eq. from 1 to m: 

Cli = ^ - YT=1 a i (/* + S i)^j = a i-l s i-l - a i (E£Ll a i (fi + S i)) 

m 

where ao = 1 — E a i an d equation on scaling parameter: 
n = nr = n (E™ i a i (fi + s i) - l ) 



Appendix C Stationary age distribution and Euler-Lotka 
equation in the continuous case 

Stationary points of this system are: 
-for age structure: 

Cil ~ lSl ~ 1 = E™ i a 'i (fi + s i) for i to m - 1 
therefore 

Elli a* (/» + s 
which implies: 

1 E™ i a * (/« + s ») 

a soSi 

a 2 = 2 

(E™ i <*i (/; + Si)) 

i-l 

a I] s j 

«i = 1 

(EiLl a i (fi + s 0) 
and ^ = 1 implies 



21 



a 



1 + 



so 



S Si 



I (this is similar to Euler-Lotka equation) 



V 



thus: 



i-l 



(E£iOi(/i + *i))' 



a 



1 + 



so 



SoSl 



i-l 

J=0 



E£i«i(/i + *0 (E£i + 



(E™ i a« (/* + s 4 )) 4 



+ ... 



therefore stable age structure is unique vector of frequencies among age 
classes, conditional on average malthusian growth rate of the population. Now 
let us proof the equivalence with Euler-Lotka equation. After substitution of 

stable age frequencies to h = n |V™. n ^— - — 1 ) equation we obtain: 

V a ) 



no = n ( £ i=0 — - 1 ) = 



1 



n 



n 



a 



V 



ao/o + 



soao/i 



soSiao/2 



V 



+ ... + 



i-l 

n Sjdo/i 

j=0 



(E™ 1 a * (/i + s i)Y 



so/i 



S0S1/2 



n sjfi 
3=0 



v 



E^iai(/i+«o (E™i^(/i+or 



(E™i a * (/< + *))' 



frequency equilibrium implies that growth rates in all age classes are equal 

to the average growth rate A = Eili a i (/» + s i) — 1, which implies: 

i-i 

, , n sjfr 

So/1 S Si/ 2 j=0 



/o + 



+ 



E^iOi(/i + *i) (YZi^Ui + si))- 

+... - 1 = ei=i (/» + s 4 ) - 1 

then 



+ ... + 



(EHi a » (/« + s *))* 



/o 



+ 



so/i 



+ 



so s 1/2 



.. + 



i-i 

n 

J=0 



(EHi a< (/i + «<)) 



i+i 
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this is Eulcr-Lotka equation 



Appendix D Impact of density dependence on the dynamics 

After substitution of the stationary age structure to the formula for turnover 
coefficient we have (note that 1 + r = EHi a i ifi + s 0) : 



E 



i=0 



i-1 

a EI s j/* 

3=0 



l-l 



E™o^n«i/i(l+r)' 

3=0 



L(f(t),s(t)) 



(1 + r) 



rn 



Em 
i=0 



(E" 1 Oi (/i + ««))' = 

i-1 ~ _ i-1 

«o n s j (! - «») Er= a o n ^ u + r ) m ~ 4 o - 



j=0 



j=0 



(E™i^ + 



(l + r) 



m 



E™on^(i+^r 

_ 3_=0 

E™ n ^-(l+rr-'a-ao 

Analogously we can substitute the stable age structure to the growth equa- 
tion (for simplicity without logistic suppresion coefficient): 



n = n 



Em 
i=l 



i-1 \ 

°0 II s 3 ifi + s i - 1) 

J'=0 



(l + r)' 



V 



) 



(l + r)' 



V 



nag 



i-1 



(l + r) 



e™ i n (i + ^ + e™ i ff (i + r) m - 1 s , - Er=i n 1 ^ (i + 

.7=0 j=0 j=0 

""" ( YZ, n (1 + r)" 1 ^ - E™1 ff Si (1 + r) m - (1 - s,: ) = 



(l + r) 



3=0 



i-1 

for simplicity denote F — E"=i II s jfi (1 + r ) m " 1 and averaged per capita 

3=0 

i-1 

mortality as T — Eti II s j (1 + r ) m * (1 ~ s *) m effect above equation have 

3=0 

form: 



n = -J---^(F-T) 



nag 



(l + r)' 



(l + r)' 



r ( - i 



(l + r) ; 



-T(L(f(t),s(t))-l) 
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F 

because L (f(t), s(t)) = — , when we remove the auxiliary notation and add 

the logistic suppression ^1 describing density dependent survival of the 

newborns, we have: 

n = na Q 1=0 (l (f(t), s(t)) (l - ^) - l) 



na (L (/(f), »(*)) (l-^)-toin *i (1 + rV (1 - *) 



when we introduce the nest site lottery by assuming that population is in 
the equilibrium size, we obtain: 



n = na [ - 1 ) YT=\ II ~ 71 (1 " «i) 

j=o(l+r) 



Lpop 



Appendix E Long age classes 

Assume that there are some number of small age classes where individuals did 
not differ in payoffs. Under this assumption we can significantly reduce number 
of age classes by substituting them by large age classes N v = n g + ... + n k 

N v = (n g + ... + n k y = n g+1 + h g+2 + n g+3 + ■■• + n m = 

k fc-1 
= J2 ( s i-l n i-l - n i) = So-1%-1 - (1 - s v ) J2 Tli - n k 

i=g i=g 

k 

= S V -lUg-l - (1 - Sy) J2 n % ~ S Vnk 

i=g 

Therefore system can be presented in the form: 



= ^N i f i -(l-s )N -s n k (59) 

i 

N v = S V -Wg-l - (1 - S V ) N v - S V Uk (60) 



Now we can assume that inside of large age class system is in equilibrium 

TL ' 

(jlip and simplify. Because a, t = — and that survivals s v are the same inside 

n 

each long age class, we have 

<ri g~ 1 

n k = ^hfv (gl) 

(E l=0 a i {fi+ S i)) 
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Now we can substitute small age class notation with long age classes. Be- 
cause:from (l61j) we have 

k -9 s i-l 

N v = n g + ... + n k = n g *^2 „ " 1 ( 62 ) 

i=l (Ei=O a i(fi + S i)) 

thus 

m 



9 k-g 

E 



i=i (J2iLo a i (ft + s dY 

which substituted to (|6"T1) leads to: 



at qk-g qk-g 

SvTlh = : — ; = i\ v - 

k—q „i— 1 /"^ m It i \\k-g k—q 

E j= " "~t a * (/i + Sl)) E sir 1 (E™ o a* (/< + si)) k - 9 ~ l 



q k-g 



i=l 

(64) 



Thus coefficient D (s v ) — describes frac- 

ff^CE^oo, (/i + *)) fc_lM 

tion of individuals passing to next large age class. In effect equation (jfJOj) can 
be presented the form: 



N v = N v -iD (s v -x)-N v (1 - s v )-N v D (s v ) = N V -%D (s v -i)-N v {l-s v + D (s v )) 

(65) 

where N v _iD (s„_x) is entering from previous age class, -/V„-D (s„) are sur- 
viving to the next age class and N v (1 — s„) die. k — g is the lenght of large age 
class and can be substituted by something. Assume that h v = k — g is the lenght 
of i>-th long age class, i.e. the number of small age classes). After simplification 
and transition into Malthusian form we obtain 

No = N f £ jifi - (1 - s ) - D ( So ) 
Nv = N v (^^D (s w _x) - (1 - s v ) - D (s v ) 



Appendix F Derivation of system a 



We start from the Malthusian sysytem describing exponential in subpopulation 
of carriers of the j-th strategy: 
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N v-1 



N l = N J [itf D ( s D - C 1 - - ^ ( s ?) J 

According to (|23p above system can be transformed into frequency replicator 
dynamics of age classes: 

4 = A T l ° (4" 1 ) + 4 M - D W - ^ti 4 (/] + 4) + D ( s ?)] 

Malthusian equation describing growth of subpopulation of J-th strategists 

= n? (e ^/i - 1 1 - *?) - ^ (*?)) +sr =1 («n - n - •?) - ^ os) 
= e Ni jfj~ N j [(! - s °j)- D («?)] +e"i at ^ (-r 1 ) -E"ii ^ (i - *?) - 

= E -^j/j - ElT=o N j i 1 - ')) + E™ i NJ-'D (sj- 1 ) - o ^ = 

i 

= E ^j/j - E™ ^ (i - - ^ (*7*) = 

= ^ f E 4/j - E™ o 4 (1 - s$ - A T D (,J 



= ^ (e 4/i + e™ o 4 s J - 4 mi? ( s r) - 1 

and now we can formulate replicator dynamics describing evolution of the 
gene pool: 



h = pj (E A )fl + EIT=o 4 S J - A T D ( S T) - E z v, [T, + E™ a;*; - a?d ( s ™) 
= » ( ( E Aj/| -E»p«E 4/i) + (E™ o 4 s ? - E, p. E™ ^) - 



therefore for system a it will have form 

4 = A 7 lD (sT 1 ) + 4 M - ^ - (E™i 4 (/i + - 



2G 



Pj = Pj ((Ij ~f) + (Sj - s) - (AfD(sf) - J2 z PzATD (s?))) 
and equation on scaling parameter 

N = N E7 =1 Pj (E™ i A) (/) + aj) - - 1) = N (/ + g - D? lo bal l) 

where fj = Eti^fj f = EUPjfj *j = E£i 44 
■5 = El"-i Pi*i ^«*oJ = Ej=i PjA™D(s™) 

Then to add the neutral density dependence the fertilities should be multi- 
plied by logistic suppression coefficient 1 — — I . 



Appendix G Derivation of system b 

System b produces more complicated form of equations. Again let us start from 
malthusian equations: 

*? = ^(E^/Ml-^ )-^ )) (66) 

NJ = NJ^-j^D(s^)-(l-s^)-D(s])j (67) 

Let us derive malthusian equations describing growth of age classes in global 
population.Then N l = E} N] and p) = j± ,since N l = E} N] 

N° = Ej N° (e ^ff (1 - S °) _ D (s o } ^ = E o ^ m _ N o {l _ s?) _ (s o^ = 

= Ej E, JVJ/f - E>° (i - -?) - E° ^ (-?) 

AT* 4» 

because £ . E* N jfi =T,i N^jPjfi = N° EiJ^EjPjfi = ^° E, ^ E „ 
above equation has form: 

^° - ^° (e, ^ E,/-;/: - E>° (i - *?) - eW («?)) 

analogously: 
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N v = Ej ^-±-D (sy 1 ) (1 - .J) - D (.J)j = 

7VTO-1 

because E^r^^r 1 ) =^- 1 E^r 1 ^( s r 1 ) =^iv^ E ^r 1 ^( s r 1 ) 

above equation has form: 

- ^ Ej pV d ( s V) Ej p? (i - »?) - Ej pP (-?)) 

therefore we obtain system: 

^° - ^(e^E^ + E^-E^^)- 1 ) ( 68 ) 

\ * 3 3 3 / 

& = ^(^E^^+E^-EpP^)- 1 ) ( M ) 

\ 3 3 3 j 

for each v 

Now we can calculate Malthusian parameter: 

M = A ° (e^E^/J + E>°*° - 1 - E>^ (,?)) + 

+ E2=i ^ E,- pj~ 1d (sir 1 ) + E, Pj-J i E,4"K)) 

= Ei ^ Y.M + A° E • P?*? - ^° - ^ E ■ P a jD (a?) + 

+ EZi a*- 1 E, pT'd (s^+YZ i ^ Ej /'>; E:", ^-E™ i ^ Ej pj^ (4) = 
= E, i' EjPjfj- TZ o ^E.^0-4)- Er=i ^ Ej p}"^ - 

e:"„-i'E,///>(^! 
- e, i' e, pj (/j + 4) - ^ m Ej - 1 

Then we can derive replicator dynamics of global population age structure: 

AV = AV Ej pVd (sr 1 ) - E,- pJ (i - «3!) - E, p]d - 
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- (E, A* V ; p« (fj + 4) - ^ V, ,,■;'!) (a?) l) 
--A*- 1 J2 j p v j - 1 D(s" j - 1 ) + 

-A v IE, P?4 - Ej P V 3 D (s)) - V, A* V, pj (/j + S j) + A m V, pJ*D 



A«- i E 3 -pr 1 -D(ar 1 ) + 



(E, p? (4 - D (*?)) - (E, ^ E, p} (/j + 4) - Am E; p^ (4?) 



then from malthusian parameters from equations (|67p and (|69[) we can derive 
dynamics of proportions of strategies inside age classes: 



'n 1 - 1 



because A i pj = NJj/N (see appendix A) 



E, 



A l p\ 

A* 
1^-1 



-D(a\ 



D^) + E j pj(i-4)+E j Pp(^)) = 



-^nm-p)?<i^h D ^ + 

E, (i - 4) • E (4) - (i - 4) - B (4) 



.4 



i-l r 



.4' 



(4 _1 ) -p1 Ei^ (4 *)H [(4 - e,p}4) - (4) - e,pP (4) 



and 



AH 



p° = p° I Ei -fijsfj - (! - 4) - ^ (4) - E. pS (E, ^ - (i - -2) - d « 



N', 



p"=p" 



(e, - (i - 4) - d (4) - E„pS (e, ^Ji - (i - -2) - d (-2))) 



p° = p" f Ei ^/j - (i - 4) - (4) - E„ E, ^rtt + E.pS (i - -S) + E.pS^ (-2) 



.4° 



P? 



A^p*- A i p i \ r / 

E! ~A0~-fj ~ P] Et, E l ~~fio~fl I +Pj (^ s i ~~ E.7 p 1 



3 ^' 3 



£(4)-E,p°£(4 
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therefore for system b: 



1 

AO 



+p) m - e:=i^4) - (d( S )) - e:=ip^(4))] 

A i = ^- 1 E^iPj£»(4" 1 ) + 

+A l 



for all j 



e;=i pH4 - D (4)) - ^ e™=i ^ (// + *j) - A m e;=i 

and equation on scaling parameter 

at = n (e™ ! ^ e;=i p)(f] + 4) - ^ m £?=i p r ) - *) 

in simplified form: 



1 

A " 



(E, - e„ e 4 ^Pi/i) + p° [(*? - 5°) - - 



i l = A i ~ 1 D i ~ 1 + A l [s l - D l - (s + f - A m D m )] 
and equation on scaling parameter 
N = N(s + f- A m D m - 1) 

where = EJ =1 p)D (s^ 1 ) & = £™ =1 P%D (*•) 

* = E™=i a v *° r = t,v=i pi ft f = e™=i A v r 



¥ - V" 1 »* s l 

° L-iv=\l J v v 



As inthe previous appendix, to add the neutral density dependence the fer- 

( N\ 

tilities should be multiplied by logistic suppression coefficient I 1 — — J . 
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Fig.l The difference between two alternative formulations of the problem: 
system a describes the evolution of the gene pool according to age structures of 
carrier subpopulations, system b describes evolution of the global age structure 
driven by strategy selection in age classes. 
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